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Introduction 

High  manoeuvrability  requirements  on  modem  helicopters  require  the  necessity  of  deeper 
understanding  of  dynamics,  and  in  turn  possibilities  for  control  of,  dynamic  stall.  Unsteady  flow 
and  resulting  aerodynamic  loads  strongly  affect  the  propulsive  efficiency  of  highly  loaded 
helicopter  rotor  blades  [Ham:  1968],  [Johnson:  1972],  Thus,  there  is  currently  much  interest  in  the 
unsteady  aerodynamics  of  flows  relevant  to  rotorcraft  (manned  and  unmanned)  air  vehicles.  The 
flowfields  in  these  applications  exhibit  unsteady  separation  followed  by  the  formation  of 
dynamic-stall-like  vortices  whose  evolution  and  interaction  with  flexible  aerodynamic  surfaces 
have  a  significant  impact  on  flight  stability  and  performance  and  in  turn  on  structural  design 
[Barwey:1994],  Analysis  of  these  flows  is  complicated  further  by  their  mixed  laminar 
transitional  turbulent  character  at  moderate  Reynolds  numbers,  as  well  as  by  the  broad  range  of 
possible  parameters,  kinematics,  and  configurations.  To  facilitate  progress  in  the  understanding 
and  prediction  of  the  relevant  fluid  dynamic  mechanisms,  it  is  natural  to  consider  methods  that 
provide  simplification  of  the  flow  phenomena  by  separating  them  into  individual  modes.  The 
technique  of  Proper  Orthogonal  Decomposition  (POD),  see  [Holmes:  1998]  is  a  popular  way  of 
accomplishing  this  task.  However,  while  POD  is  capable  of  extracting  the  most  energetic  parts  of 
the  flow  field,  it  has  been  shown  to  lack  ability  of  highlighting  subtle,  oscillatory  phenomena  that 
are  nevertheless  implicated  in  important  physical  processes  such  as  the  shear  layer  separation. 

Unsteady  flows  over  plunging  and  pitching  airfoils  with  large  excursions  in  effective  angle  of 
attack  exhibit  the  phenomenon  termed  dynamic  stall,  a  process  characterized 
phenomenologically  by  unsteady  separation  and  by  the  formation  of  large-scale  leading-edge  and 
trailing-edge  vortices,  which  exert  difficult-to-predict  variations  in  aerodynamic  loads.  Distinct 
stages  of  the  evolution  of  dynamic  stall  are  a)  vortex  formation,  b)  vortex  convection,  c)  stall 
onset  and  d)  stalled  stage.  Comprehensive  reviews  of  this  phenomenon,  first  discovered  and 
studied  extensively  in  the  context  of  helicopter  rotor  blades,  has  been  given  by  McCroskey 
[Mccroskey:1982],  [Carr:2012],  and  [Ekaterinaris:1998]. 
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Aimdyn  and  AFRL  have  been  studying  physical  mechanisms  of  dynamic  stall  using  the  methods 
of  Koopman  operator  theory.  The  initial  seed  research  has  suggested  existence  of  high-frequency 
effects  in  leading  edge  vortex  shedding  events,  not  present  in  the  static  stall  case.  The  research 
guided  by  dynamical  systems  methods  and  efficient  numerical  methods  for  Koopman  Mode 
Decomposition  already  lead  to  interesting  outcomes.  Among  those  are  insight  into  the  physical 
effects  of  the  oscillatory  nature  of  the  incoming  velocity  that  the  wing  experiences  due  to  pitch  or 
plunging.  For  example,  the  examination  of  a  model  system  -  a  cylinder  in  an  incoming 
oscillatory  flow  -  showed  some  of  the  physical  effects  (e.g.  broadening  of  the  spectrum) 
observed  in  the  pitching  airfoil  case  are  present  in  the  model  system  and  shed  light  on  the 
dynamics  of  aspects  of  pitching  airfoil  dynamics  as  a  nonlinear  interaction  of  the  forcing  by  the 
oscillating  incoming  flow  frequency  and  natural  frequency  of  vortex  shedding  dynamics.  In 
collaboration  with  the  Army  Research  Laboratory's  Dr.  Bryan  Glaz,  we  found  the  precise  state- 
space  description  of  the  mechanism  of  the  broadening  of  the  spectrum  by  examining 
perturbations  of  the  limit  cycling  motion  by  external  forcing  frequency. 


Koopman  Modes  for  Stationary  Airfoil  with  Different  Degrees  of  Oscillation  Amplitude 

Koopman  mode  decomposition  is  based  on  the  surprising  fact,  discovered  in  [Mezic:  2005],  that 
normal  modes  of  linear  oscillations  have  its  natural  analogue  -  Koopman  modes  -  in  the  context 
of  nonlinear  dynamics.  To  pursue  this  analogy,  one  must  change  the  representation  of  the  system 
from  the  state-space  representation  to  the  dynamics  governed  by  the  linear  Koopman  operator 
([Koopman:  1931])  on  an  infinite-dimensional  space  of  observables.  Contrary  to  the  proper 
orthogonal  decomposition,  the  dynamic  mode  decomposition  contains  not  only  information 
about  coherent  structures,  but  also  about  their  temporal  evolution. 


Based  on  snapshots  of  the  flow,  we  can  approximate  the  Koopman  Modes  using  an  Amoldi-like 
algorithm  sometimes  called  dynamic  mode  decomposition  (DMD)  ([Schmid:  2008],  [Rowley: 
2009])  which  computes  eigenvalues  based  on  the  so-called  companion  matrix. 


Given  a  sequence  of  equispaced  in  time  snapshots  from  numerical  simulations  or  physical 
experiments,  with  At  being  the  time  interval  between  snapshots,  a  data  matrix  is  formed  with 
columns  that  represent  the  individual  data  samples  Uj  6  Rn,j  =  0, ... ,  m,  with  j  representing  time 
jht.  The  companion  matrix  is  then  defined  as: 


C  = 


where  q,  i  =  0, ... ,  m  —  1  are  such  that: 


(  0  0  •  •  •  0  co 

10  0  Cl 

0  1  0  c2 

V  o  0  •  •  •  1  cm— i 


m—  1 

um  =)  CiUj  +  r 
3=  0 


\ 


/ 


and  r  is  the  residual  vector. 


The  spectrum  of  the  Koopman  operator  restricted  to  the  subspace  spanned  by  Uj  is  equal  to  the 
spectrum  of  the  infinite-dimensional  companion  matrix  and  the  associated  Koopman  modes  are 
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given  by  Ka  (provided  that  a  does  not  belong  to  the  null  space  of  K),  where  K  = 
[u0,Wi,  is  the  column  matrix  (vector-valued)  of  observables  snapshots  at  times 

0,  At, (m  —  l)At  and  a  is  an  eigenvector  of  the  shift  operator  restricted  to  Krylov  subspace 
spanned  by  which  the  Companion  matrix  is  an  approximation  of.  The  approximate  Koopman 
eigenvalues  and  eigenvectors  obtained  by  the  Amoldi's  algorithm  are  called  Ritz  eigenvalues  and 
eigenvectors. 

The  standard  Amoldi-type  algorithm  to  calculate  the  Ritz  eigenvalues  Xj  and  eigenfunctions  Vj  is 
as  follows: 


1. 

2. 


Define  if  =  [u0,u1( 


Find  constants  Cj  such  that: 


771—1 


r  =  Um  CjUj  =  um  ~  Kc,  r±span{u0, ...,  um_i}. 

i= 0 

This  can  be  done  by  defining  c  =  K+um,  where  K+  is  the  pseudo  inverse  of  K. 

3.  Define  the  companion  matrix  C  as  above  and  find  its  eigenvalues  and  eigenvectors: 

C  =  T-'ACA  =  diag(\u--  ■  ,Xm), 

where  eigenvectors  are  columns  off-1.  Note  that  the  Vandermonde  matrix  f 


T  = 


1  Aj  A? 

1  A2  \\ 


Ai  \ 

\  771—1 


1  Am  A^  •••  A--1/ 

diagonalizes  the  companion  matrix  C  ,  as  long  as  the  eigenvalues  Xlt ... ,  Xm  are  distinct. 


4. 


Define  Vj  to  be  the  columns  of  V  =  KT  1. 


Then,  the  Amoldi-type  Koopman  Mode  Decomposition  gives: 

771 

Vfc  =  [0, 1,  •  •  •  ,  m],  uk  =  ^2  (:-  3) 

i= i 


We  studied  the  stationary  airfoil  simulation  provided  by  our  collaborator  Dr.  Bryan  Glaz  at  the 
U.S.  Army  Research  Laboratory.  The  simulation  was  for  a  statically  pitched  NACA  0015  airfoil. 
The  airfoil  chord  was  0.135  m,  fresstream  velocity  magnitude  was  100  m/s,  and  the  freestream 
was  oriented  at  a  static  18  degree  angle  of  attack.  We  also  studied  the  case  with  6  deg  oscillation 
amplitude  as  well  as  2  deg,  4  deg  and  8  deg  oscillation  amplitude. 

Figure  1  shows  Fourier  and  KMD  spectrum  for  u-velocity. 
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Figure  1.  Fourier  and  KMD  spectrum.  FFT  is  in  black;  KM  spectrum  is  red. 


The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  2,  DMD  spectrum  in  frequency 
w,  exponential  rate  mu  plane  is  shown  in  Figure  3. 
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Figure  2.  Koopman  Mode  Eigenvalues. 
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Figure  3.  DMD  spectrum  in  frequency  w, 
exponential  rate  mu  plane. 


The  recomposition  of  signal  using  40  first  pairs  of  the  highest  magnitude  modes  sorted  by  abs  of 
norm  of  Vj  has  been  performed.  Figure  4  shows  the  reconstructed  signal  (red)  vs  the  original 
data  (blue). 


Figure  4.  The  reconstructed  signal  (red)  vs  the  original  data  (blue)  for  the  random  location. 


Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  mode  corresponding  to 
frequency  470  Hz  are  shown  in  Figure  5. 
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Figure  5.  Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  mode  corresponding 

to  frequency  470  Hz. 


Figure  6  shows  Fourier  and  KMD  spectrum  for  u-velocity  for  the  airfoil  case  with  6  deg 
oscillation  amplitude. 
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Figure  6.  Fourier  and  KMD  spectrum  for  the  airfoil  case  with  6  deg  oscillation  amplitude.  FFT  is 

in  black;  KM  spectrum  is  red. 


Figure  7  shows  Fourier  and  KMD  spectrum  for  u-velocity  for  the  airfoil  case  with  2  deg 
oscillation  amplitude. 
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frequency  wj 

Figure  7.  Fourier  and  KMD  spectrum  for  the  airfoil  case  with  2  deg  oscillation  amplitude.  FFT  is 

in  black;  KM  spectrum  is  red. 

Figure  8  shows  Fourier  and  KMD  spectrum  for  u-velocity  for  the  airfoil  case  with  4  deg 


in  black;  KM  spectrum  is  red. 

The  real  part  for  the  Koopman  mode  corresponding  to  frequency  24  Hz  for  the  u-velocity  for  the 
airfoil  case  with  6  deg  oscillation  amplitude  is  shown  in  Figure  9 A,  the  real  part  for  the 
Koopman  mode  corresponding  to  frequency  23.5  Hz  for  the  u-velocity  for  the  airfoil  case  with  2 
deg  oscillation  amplitude  is  shown  in  Figure  9B,  the  real  part  for  the  Koopman  mode 
corresponding  to  frequency  23.5  Hz  for  the  u-velocity  for  the  airfoil  case  with  4  deg  oscillation 
amplitude  is  shown  in  Figure  9C. 
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Real  Part:  Freq=23.5661 


0  0.1  0.2 


Figure  9 A.  The  real  part  for 
the  Koopman  mode 
corresponding  to  frequency 
23.5  Hz  for  the  airfoil  case 
with  6  deg  oscillation 
amplitude. 


Real  Part:  Freq=23.5208 


0  0.1  0.2 


Figure  9B.  The  real  part  for 
the  Koopman  mode 
corresponding  to  frequency 
23.5  Hz  for  the  airfoil  case 
with  2  deg  oscillation 
amplitude. 


Real  Part:  Freq=23.5955 


Figure  9C.  The  real  part  for 
the  Koopman  mode 
corresponding  to  frequency 
23.5  Hz  for  the  airfoil  case 
with  4  deg  oscillation 
amplitude 


Connection  between  a  basic  version  of  the  dynamic  mode  decomposition  and  the  discrete¬ 
time  version  of  the  generalized  Laplace  analysis 

We  studied  the  connection  between  a  basic  version  of  the  dynamic  mode  decomposition  and  the 
discrete-time  version  of  the  generalized  Laplace  analysis  to  apply  to  the  airfoil  case. 

Theorem:  Let  f(x;z)  be  a  field  of  observables  f(x;z) :  M  x  A  ->  Rm,  where  the  observables  are 
indexed  over  the  set  A.  We  will  occasionally  drop  the  dependence  on  the  state-space  variable  x 
and  denote  f(x;z)  =  f(z)  and  the  iterates  of  f  by  f(T‘x;z)  =  f(z).  Let  Ao, ...,  A  be  the  eigenvalues  of 
Uc  such  that  M>IMI>  ...  >  INI-  Then,  the  Koopman  mode  associated  with  Ak  is  obtained  by 
computing 


We  showed  an  explicit  relationship  between  a  basic  version  of  the  Dynamic  Mode 
Decomposition  (DMD)  and  the  Koopman  Mode  Decomposition  (KMD)  of  dynamical  systems, 
that  allows  for  estimates  of  validity  of  approximation  of  Koopman  modes  by  DMD  modes.  We 
also  linked  the  recently  introduced  Generalized  Laplace  Analysis  and  the  inverse  of  the 
Vandermonde  matrix.  We  proved  the  equivalence  between  dynamic  modes  and  Koopman  modes 
over  infinite  iteration  of  the  dynamical  system. 

Prony  Methods 

We  also  studied  Prony  methods  [Plonka:  2014]  for  recovery  of  structured  functions  and 
implemented  the  following  algorithm  in  C++. 
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Algorithm  (Classical  Prony  method) 


Input: 

McN,  sampled  values  h(k),  k  =  0,. .  .,2M-1,  of  the  exponential  sum 

M 

h(x )  :=  cj  e^jX  i  x  >  0 , 
i= i 

1 .  Solve  the  following  linear  system 


where 


H„(0)  (w)";1 


—  {h(M  +  m)) 


M- 1 

m=0 


Hm(0) 


Mo) 

Mi) 


Mi) 

M2) 


h(M  —  1)  \ 
h(M) 


+  m)) 


M  —  1 
k^m— 0 


h(M-  1)  h(M) 


h(2M  —  2)  / 


2.  Compute  all  zeros  Zj  c  D,  j  =  1,. .  ,,M,  of  the  Prony  polynomial 

M  Af— 1 

p{z )  :=  -  Zj)  =  ^  Pfc  +  zM,  zee 

j=l  fc=0 

i.e.,  calculate  all  eigenvalues  of  the  associated  companion  matrix 


(  0 

0  . , 

..  0 

Po 

1 

0  . , 

..  0 

Pi 

c M(p)  := 

0 

1  . 

..  0 

-P2 

l  0 

0  . . 

. .  1 

—PM-1 

and  form  fj  =  log  Zj  for  j  =  1,. .  ,,M,  where  log  is  the  principal  value  of  the  complex  logarithm. 

3.  Solve  the  Vandermonde  system 

Vm(z)  (cjOJC  -  . 

Output: 

fj  £  [ — Ct ,  0]  +  i  [ — 7T,  7r),  Cj  £  C,  j  —  1.  .  .  .  ,  Al. 


Bifurcation  around  a  Critical  Reynold’s  Number 

We  were  looking  at  the  theoretical  explanation  of  our  numeric  observations  of  bifurcation  around 
a  critical  Reynold’s  number.  Reading  paper  [Bagheri:  2013],  we  realized  that  we  could  look  at 
this  issue  as  a  forced  oscillator  set  up  similar  to  the  one  in  paper  [Vance:  1989],  We  were 
studying  this  paper  to  apply  the  theory  to  our  case. 
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We  compared  analytic  solutions  for  the  bifurcation  sets  of  the  truncated  normal  form  of  the 
dynamic  equations  for  p=l  (Figure  10 A),  p=2  (Figure  10B)  and  p=40  (Figure  10C). 


Figure  10 A.  p=l 


Figure  10B.  p=2 
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Figure  IOC.  p=40 

We  researched  the  following  papers:  [Wang:  2003],  [Lin:  2008],  The  main  idea  there  is  that  for 
large  periods  T  of  perturbation  chaotic  behavior  arises  in  systems  with  Hopf  bifurcation.  This  is  a 
restricted  result  in  that  the  considered  perturbation  consists  of  kicks.  We  are  considering 
pursuing  simulations  of  a  simple  system  where  the  perturbation  is  smooth,  and  pursuing  CFD 
with  oscillating  forcing  in  form  of  kicks  with  very  large  period. 

The  team  members  studied  the  cylinder  wake  simulation  provided  by  our  collaborator  Dr.  Bryan 
Glaz  at  the  U.S.  Army  Research  Laboratory.  There  are  two  cases  for  p=l  and  for  p=2  (assuming 
q=2  in  the  paper  [Vance:  1989]).  For  p=l  the  cylinder  is  oscillated  at  twice  the  natural 
bifurcation  frequency,  and  at  p=2  it  is  oscillated  at  the  natural  shedding  frequency  of  -39  Hz. 
Within  each  p,  there  are  large  amplitude  oscillation  and  small  amplitude  oscillation.  For  p=2,  the 
tecplot  output  files  were  output  at  every  12  time  steps  (CFD  time  step  size  =  le-4  s).  The  first 
5000  time  steps  and  the  corresponding  output  files  are  for  the  stationary  cylinder.  Then,  the 
cylinder  is  oscillated  for  the  remaining  time  steps.  After  5000  time  steps,  each  case  is  run  enough 
time  to  capture  about  10-12  periods  of  the  oscillation  frequency.  Everything  the  same  is  for  p=l, 
except  files  are  outputted  at  every  6  time  steps.  It  was  done  because  the  cylinder  oscillates  twice 
as  fast  as  the  p=2  case  and  to  make  sure  to  capture  higher  frequency  content. 

Figure  1 1  shows  Fourier  and  KMD  spectrum  for  u- velocity  for  p=l  and  small  amplitude. 
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Figure  1 1 .  Fourier  and  KMD  spectrum.  FFT  is  in  black;  KM  spectrum  is  red. 


The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  12,  DMD  spectrum  in  frequency 
w,  exponential  rate  mu  plane  is  shown  in  Figure  13. 


tigenvaiues  on  tne  unit  uircie  (Normalized) 


witnout  mean  removal,  normalized 

12 
10 
8 
6 
4 
2 

'"-Tooo  -500  o  500 

Figure  13.  DMD  spectrum  in  frequency  w, 
exponential  rate  mu  plane. 


The  recomposition  of  signal  using  20  first  pairs  of  the  highest  magnitude  modes  sorted  by  abs  of 
norm  of  Vj  has  been  performed.  Figure  14  shows  the  reconstructed  signal  (red)  vs  the  original 
data  (blue). 

0.56 
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14.  The  reconstructed  signal  (red)  vs  the  original  data  (blue)  for  the  random  location. 

shows  Fourier  and  KMD  spectrum  for  u-velocity  for  p=2  and  small  amplitude  by  using 
not-normalized  Vandermonde  matrix  (top)  and  normalized  Vandermonde  matrix  (bottom). 


Figure 
Figure  15 
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Figure  15.  Fourier  and  KMD  spectrum.  FFT  is  in  black;  KM  spectrum  is  red. 


The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  16,  DMD  spectrum  in  frequency 
w,  exponential  rate  mu  plane  is  shown  in  Figure  17. 
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Figure  17.  DMD  spectrum  in  frequency  w, 
exponential  rate  mu  plane. 


The  recomposition  of  signal  using  20  first  pairs  of  the  highest  magnitude  modes  sorted  by  abs  of 
norm  of  Vj  has  been  performed.  Figure  18  shows  the  reconstructed  signal  (red)  vs  the  original 
data  (blue). 


Figure  18.  The  reconstructed  signal  (red)  vs  the  original  data  (blue)  for  the  random  location. 

Reduced  Order  Model  (ROM)  of  Dynamic  Stall 
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The  team  members  began  working  on  the  Reduced  Order  Model  (ROM)  of  Dynamic  Stall  using 
skew-projection  methods  used  in  conjunction  with  the  Koopman  Mode  Decomposition. 

We  used  the  following  four  cylinder  simulations  provided  by  Dr.  Bryan  Glaz  at  the  U.S.  Army 
Research  Laboratory.  The  cylinder  was  0.002m  in  diameter.  Each  simulation  was  run  for  6  s. 

The  number  of  nodes  in  the  grid  is  29954.  For  the  stationary  cylinder  case,  there  are  output  files 
at  every  12  time  steps,  so  the  delta  t  between  outputs  is  0.0012  s.  For  the  oscillating  cylinder 
cases,  flow  snapshots  were  outputted  starting  at  the  25th  time  step,  up  to  the  60,000th  time  step 
in  increments  of  25  time  steps  (each  time  step  is  le-4  s).  From  0.0s  to  0.5s  (i.e.  5000  iterations) 
Reynolds  number  (Re)  was  58.3  that  corresponds  to  incoming  velocity  of  0.5  m/s.  The  critical  Re 
for  a  cylinder  is  about  Re~40  when  it  starts  exhibiting  the  von  Karman  wake  instability.  So  the 
interval  from  0  to  0.5  s  at  Re=58.3  corresponds  to  the  Landau  equation  when  the  external  input  is 
greater  than  the  bifurcation  value.  Then,  starting  at  0.5  seconds,  the  Re  is  being  oscillated  by 
oscillating  incoming  velocity.  The  oscillating  Re  corresponded  to  Re  =  58.3  + 
35*sin(2pi*omega*t),  where  omega  corresponded  to  2  Hz.  The  reason  35  was  selected  for  the 
oscillating  Re  amplitude  is  because  a  large  enough  amplitude  was  needed  such  that  the  Re  would 
oscillate  above  and  below  the  critical  value  predicted  by  the  regular  Landau  equation.  The  reason 
2  Hz  was  selected  is  that  it's  an  order  of  magnitude  slower  than  the  von  Karman  vortex 
frequency,  which  is  about  39  Hz  for  Re  =  58.3.  Two  simulations  for  the  oscillating  cylinder  case 
were  performed  for  the  case  of  Re  =  58.3  +  5.85*sin(2pi*omega*t)  and  Re  =  58.3  + 

1  *sin(2pi*omega*t). 

We  used  the  following  algorithm: 

Let  Xj,  j=l,2,...,  4969,  be  eigenvalues  corresponding  to  Koopman  Modes  Vj,i,  where 
j=l,2,...,4969,  coordinates  i=l,2,...,29954,  for  the  stationary  cylinder  case.  We  find  the  Koopman 
Mode  V ind,i  from  the  stationary  cylinder  case  that  corresponds  to  the  frequency  of  39.26  Hz 
(Figure  19). 


0.15 


o  Freq=39.26  Hz 


#7*®  Fourier 

O  Koopman  using  normalized  Vandermonde 
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frequency  wj 

Figure  19.  Fourier  and  KMD  spectrum  for  u-velocity.  FFT  is  in  black;  KM  spectrum  is  red. 

Then  we  find  adjoint  matrix  B=(inv(V*V')*V)',  that  satisfies  the  condition  that  rBVlkk=l  and 
[BV]k,i=0  for  k  ±  1. 


Let  X  of  size  4969x29954  be  the  u-velocity  corresponding  to  the  stationary  cylinder. 


13 


Contract  Number:  W91 1NF-14-C-0102 
Proposal  Number:  P-65502-EG 


Then  we  calculate  Cind,k  in  the  following  way:  for  each  k=l:size(B,l)  for  each  i=l:size(X,l)  we 
calculate  Cind,k(i,:)=dot(X(i,:),B(kv))*Vind(:).  Then  we  find  norm(abs(cind,k))  (Figure  20). 
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Figure  20.  Norm  of  Cind,k. 
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In  Figure  22  we  present  Cind,ind  vs  simulation  time  and  its  embedding  in  time  for  w=39.26  Hz  for 
Point  5  (Figure  21). 


Figure  21.  Points  around  the  cylinder. 


Figure  22.  Cind,ind  vs  simulation  time  (left)  and  its  embedding  in  time  (right)  for  Point  5  for 

w=39.26  Hz. 
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We  also  calculate  Cind,k  for  the  oscillating  cylinder  cases  by  using  the  same  algorithm  as  above 
with  the  difference  that  X  of  size  2400x29954  is  the  u-velocity  corresponding  to  the  oscillating 
cylinder  case.  In  this  case  norm(abs(cind,k))  is  shown  in  Figures  23,  24  and  25. 
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Figure  23.  Norm  of  Cind,k  for  the  oscillating  cylinder  case  with  Re=35. 
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Figure  24.  Norm  of  Cind,k  for  the  oscillating  cylinder  case  with  Re=5.83. 
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Figure  25.  Norm  of  Cmd,k  for  the  oscillating  cylinder  case  with  Re=l. 


In  Figures  26,  27  and  28  we  present  Cind,ind  vs  simulation  time  and  its  embedding  in  time  for  Point 
5  (Figure  21)  for  w=39.26  Hz  for  all  oscillating  cylinder  cases. 
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Figure  26.  Cind,ind  vs  simulation  time  (left)  and  its  embedding  in  time  (right)  for  Point  5  for 
w=39.26  Hz  for  the  oscillating  cylinder  case  with  Re=35. 


Figure  27.  Cind,ind  vs  simulation  time  (left)  and  its  embedding  in  time  (right)  for  Point  5  for 
w=39.26  Hz  for  the  oscillating  cylinder  case  with  Re=5.83. 
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Figure  28.  Cind,ind  vs  simulation  time  (left)  and  its  embedding  in  time  (right)  for  Point  5  for 
w=39.26  Hz  for  the  oscillating  cylinder  case  with  Re=l. 

Koopman  mode  projections  of  cylinder  in  oscillating  inlet  flow  show  an  interesting  transition 
that  happens  with  the  increase  of  Reynolds  number.  Specifically,  the  spectrum  is  strongly  peaked 
at  low  Reynolds  number,  and  shows  continuous  elements  as  Reynolds  number  increases.  The 
transition  starts  with  broadening  of  the  spectrum  around  harmonics  of  vortex  shedding  frequency 
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and  leads  to  merging  of  spectrum  into  a  continuous-looking  spectrum  at  higher  Reynolds 
numbers.  We  researched  connections  of  these  transitions  to  normal  form  equations  that  extend 
those  obtained  for  Hopf  bifurcation.  We  found  that  adding  a  constant  forcing  term  to  the  Hopf 
bifurcation  normal  form,  along  the  lines  suggested  by  the  paper  of  [Tsarouhas:  1987]  does  not 
lead  to  a  good  match  with  observed  transition  in  the  fluid-mechanical  problem.  Further 
investigation  lead  to  suggestion  of  adding  a  term  proportional  to  the  radial  coordinate  in  the 
normal  form,  similar  to  the  suggestion  made  in  papers  of  Young  and  collaborators  (2003,  2008) 
This  approach  proved  very  promising,  with  spectral  results  matching  closely  the  progression  of 
spectral  changes  in  the  flow,  especially  around  shedding  frequency.  The  normal  form  model, 
however,  does  not  match  the  spectrum  at  the  inlet  flow  oscillating  frequency.  Our  conjecture  is 
that  this  might  be  due  to  absence  of  pure  (radial  coordinate  independent)  oscillating  term  in  the 
current  version  of  the  normal  form. 


Conclusions 

We  believe  that  the  given  study  can  be  extended  to  capture  aspects  of  a  number  of  different 
applied  problems  of  Army  interest  in  flow-structure  interaction,  where  the  flow  affecting  the 
structure  is  unsteady  or  the  structure  is  moving  in  an  unsteady  manner. 
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